
Job: betamc.std.prg (/home/gov/faculty/ppaolino/parr)
Sun Jun 17 01:01:38 2001, pid=18242, uid=1141(ppaolino)

cheng alg used if cheng=1, johnk otherwise, cheng=       0.0000000 
parm1 =      -0.80000000       0.30000000 
parm2 =      -0.40000000      -0.30000000 
p =       0.47521835 
q =       0.69638213 
true first difference (mean)=      0.28784477 
true first difference (variance)=     0.029104255 

descriptive stats for standard beta parameters

     -0.79297857      0.058045483 
      0.29847112      0.054676967 
     -0.39265364      0.063636594 
     -0.30195275      0.061685258 
 
standard beta mse=      0.10343400 
 
alternative beta mse=      0.10349421 
 
ols mse=      0.10354559 
 
logit trans ols mse=      0.12260043 
 
het ols mse=      0.10355423 
 
logit trans het ols mse=      0.12210827 
 
log trans het ols mse=      0.13696488 
 
first differences for mean and variance for std beta (sim)
      0.28714731      0.027252838 
     0.028813027     0.0079650614 
 
first differences for mean and variance for std beta (nosim)
      0.28757435      0.027272080 
     0.029019468     0.0080051867 
 
first differences for mean and variance for alt beta (sim)
      0.27230954      0.024197912 
     0.027180189     0.0077730299 
 
first differences for mean and variance for alt beta (nosim)
      0.27267089      0.024216149 
     0.027369269     0.0077918663 
 
first differences for mean ols
      0.27220545      0.026175003 
 
first differences for mean ols (logit trans)
      0.55591324      0.052470847 
 
first differences for mean and variance for het normal
      0.26869370      0.025271372 
     0.027754338      0.012264690 
 
first differences for mean and variance for logit trans.
      0.54943712      0.052556048 
      0.39509211      0.078818106 
 
first differences for mean ols (log trans)
      0.25410209      0.028116084 
 
efficiency of standard beta vs. alternative beta=       104.79597 
 
efficiency of beta vs. ols=       111.81329 
 
efficiency of beta vs. ols (log trans)=       1002.0271 
 
efficiency of standard beta vs. het normal=       116.28142 
 
efficiency of standard beta vs. het normal (log trans)=       978.78412 
 
efficiency of standard beta vs. OLS (log trans)=       161.08833 
 
efficiency of standard beta vs. alternative beta=       99.715745 
 
efficiency of standard beta vs. het normal=       154.12679 
 
efficiency of standard beta vs. het normal (log trans)=       4678.6753 
 
standard beta parms' overconfidence=
       104.97806 
       99.171521 
       102.12402 
       101.59049 
 
alternative beta parms' overconfidence=
       99.620337 
       94.065693 
       103.28134 
       96.658827 
 
ols parms' overconfidence=
       100.05648 
       90.626171 
 
ols (log. trans) parms' overconfidence=
       100.13480 
       116.87610 
 
ols (log trans) parms' overconfidence=
       101.80935 
       92.840116 
 
heteroskedastic normal parms' overconfidence=
       100.88836 
       88.010416 
       71.838055 
       85.645750 
 
heteroskedastic log. trans. parms' overconfidence=
       99.253535 
       115.86594 
       163.36912 
       184.12290 
 
standard beta's mean first difference overconfidence
       100.56276 
 
standard beta's variance first difference overconfidence
       100.49215 
 
alternative beta's mean first difference overconfidence
       94.246106 
 
alternative beta's variance first difference overconfidence
       97.203684 
 
check standard beta's percentage within 95% conf int

      0.94100000 
      0.95300000 
      0.94700000 
      0.94700000 
 
check alternative beta's percentage within 95% conf int

      0.94900000 
      0.96300000 
      0.94300000 
      0.95200000 
 
check ols percentage within 95% conf int

      0.94700000 
      0.96500000 
 
check logit trans/ols percentage within 95% conf int

      0.93600000 
      0.90500000 
 
check het ols percentage within 95% conf int

      0.94500000 
      0.96500000 
      0.99400000 
      0.97800000 
 
check logit trans/ols percentage within 95% conf int

      0.94000000 
      0.90800000 
      0.75100000 
      0.70900000 
 
check logit trans/ols percentage within 95% conf int

      0.94000000 
      0.96300000 
 
check percentage of sig parms (standard beta)
       1.0000000 
      0.95100000 
 
check percentage of sig parms (alternate beta)

       1.0000000 
       1.0000000 
      0.80900000 
      0.14000000 
 
check percentage of sig parms (basic ols)

       1.0000000 
       1.0000000 
 
check percentage of sig parms (logit trans ols)

       1.0000000 
       1.0000000 
 
check percentage of sig parms (het ols)

       1.0000000 
       1.0000000 
       1.0000000 
      0.45600000 
 
check percentage of sig parms (het ols - logit trans)

       1.0000000 
       1.0000000 
       1.0000000 
      0.61400000 
 
check percentage of sig parms (ols -- log tranform)

       1.0000000 
       1.0000000 
 
descriptive stats for y's for successful runs

      0.41310658 
      0.34986845 
   3.0936427e-07 
      0.99997013 
descriptive stats for y's for failed runs
       0.0000000 
descriptive stats for x's for successful runs

     0.029825793 
       1.0227222 
      -3.2389446 
       3.4797200 
descriptive stats for x's for failed runs
       0.0000000 
number of failed runs
       0.0000000 

code used to generate the above results:

library cml;
#include cml.ext;
CMLset;

fails=0;
crash=0;
ystatf=0;
xstatf=0;

/* load x matrix based upon number of observations */

x = ones(500,1)~rndn(500,1);
n = rows(x);

/* produce matrices for simulated first differences */

meanx=meanc(x); @ create vector of means for the indep vars @
stdx=stdc(x);   @ create vector of stdev for the indep vars @
minx=minc(x);
maxx=maxc(x);

/* start obtaining the first difference vectors (means) */

sw=zeros(1,cols(x)-1)|eye(cols(x)-1);
xl=meanx-stdx.*sw;           @ these are kXk-1 vectors @
xh=meanx+stdx.*sw;       
xmin=meanx.*(1-sw)+minx.*sw;  
xmax=meanx.*(1-sw)+maxx.*sw;

/* start creating Monte Carlo data */
k = cols(x);
parm1 = {-0.8, 0.3};
parm2 = {-0.4, -0.3};
cheng=0;
print "cheng alg used if cheng=1, johnk otherwise, cheng=";;cheng; 
p = exp(x*parm1);
q = exp(x*parm2);

/* determine true first differences */

/* get  alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alt=exp(xl'*parm1);                          @ nxk-1 matrix @
         blt=exp(xl'*parm2);
         aht=exp(xh'*parm1);
         bht=exp(xh'*parm2);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlt=alt./(alt+blt);                    @ nxk-1 matrix @
         meanht=aht./(aht+bht);
         varlt=(alt.*blt)./(((alt+blt)^2).*(alt+blt+1));
         varht=(aht.*bht)./(((aht+bht)^2).*(aht+bht+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmt=meanht-meanlt;
         fdvt=varht-varlt;

/* print data parameters for simulation */

print "parm1 = ";; parm1';
print "parm2 = ";; parm2';
print "p = ";; meanc(exp(x*parm1));
print "q = ";; meanc(exp(x*parm2));
print "true first difference (mean)=";; fdmt;
print "true first difference (variance)=";; fdvt;

/* start Monte Carlo trials */
trial=1;
do while trial<=1000; 

print "MC trial=";; trial;

   if trial==1;
      ret=0;
      if fails==0;
         bt=0; 
      endif;
   endif;

   if fails==0;
   success=0;

/* comment out random number generator if p,q<1 */
if cheng==1;
/* create J&K (Chen) alpha */
ralpha=p+q;

/* create J&K (Chen) beta */
p1=p.>1;
q1=q.>1;
rbeta=real((p1.*q1).*(sqrt((ralpha-2)./(2.*p.*q-ralpha))) +
      (1-p1.*q1).*(1/(minc((p'|q')))));

/* create J&K (Chen) gamma */
rgamma=p+1/rbeta;

/* step 1 in J&K (Chen) RNG */
u1=rndu(n,1);
u2=rndu(n,1);
v=rbeta.*ln(u1./(1-u1));
w=p.*exp(v);

/* step 2 in J&K (Chen) RNG */

done=0;
      do while done<n;
         d=(ralpha.*ln(ralpha./(q+w))+rgamma.*v-1.3862944).<ln((u1^2).*u2); 
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v=(1-d).*v + d.*rbeta.*ln(u1./(1-u1));
         w=(1-d).*w + d.*p.*exp(v);
         done=n-sumc(d);
      endo;

/* step 3 in J&K (Chen) RNG */
y=w./(q+w); 
endif;

/* discretize y 

y1=yc.>0;
y2=yc.>0.14;
y3=yc.>0.28;
y4=yc.>0.42;
y5=yc.>0.56;
y6=yc.>0.70;
y7=yc.>0.84;

y=(y1+y2+y3+y4+y5+y6+y7)/8; */

/* comment out random number generator if p or q>1 */
if cheng==0;
/* step 1 in Johnk's RNG */

u1=rndu(n,1);
u2=rndu(n,1);

v1=u1^(1/p);
v2=u2^(1/q);

/* step 2 in Johnk's RNG */

w=v1+v2;

/* step 3 in Johnk's RNG */

done=0;
      do while done<n;
         d=w.>1;
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v1=(1-d).*v1 + d.*u1^(1/p);
         v2=(1-d).*v2 + d.*u2^(1/q);
         w=v1+v2;
         done=n-sumc(d);
      endo;

/* step 4 in Johnk's RNG */
y=v1./w;
endif;

/* cheat function to re-set Ys that equal 1 --- would choke llik function
   otherwise */
d=y.==1;
print "cheat=";;sumc(d);
y=(1-d).*y + d.*.99999;  


/* return the monte carlo data set */
    data = y~x[.,2:cols(x)];

   endif;

   if fails==0;

      do while success==0 and fails<=30;

/* maximum likelihood function for estimating parameters alpha and beta */
__title=" ** A Beta Model, Analytic Solution (Standard Estimation) **";

stval=1.0|ones(k-1,1)*0.5|1.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

       if fails>=1 and fails<=40;
          stval=(bt')/2;
          print "retrying run";
       endif;

print "mean y, max y, min y, std y";
meanc(y)~maxc(y)~minc(y)~stdc(y);

proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

proc gradproc(theta,ds);
    local y,x,z,be,h,xb,p,q,per1,gr1,gr2,grad;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    per1=p+q;
    gr1= psi(p+q).*x.*p-psi(p).*x.*p+x.*p.*ln(y);
    gr2= psi(p+q).*x.*q-psi(q).*x.*q+x.*q.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);
    endp;

proc loglik(theta,ds);
    local y,x,z,be,h,xb,p,q,l1,l2,l3,l4,l5,llik;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradproc;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400;

   if fails.>=5 and fails.<=10;
     _cml_DirTol=.000001;
   endif;
   if fails.>10 and fails.<=15;
     _cml_DirTol=.0000001;
   endif;
   if fails.>15 and fails.<=20;
     _cml_DirTol=.00000001;
   endif;
   if fails.>20;
     _cml_DirTol=.00001;
   endif;

    {b,logl,g,vc,ret}=CMLprt(CML(data,0,&loglik,stval));

/* check that parameters are reasonable */
   bt=(b');

/* determine success or failure of the trial */
   cor=vc[1,1]/vc[1,1];
         if ret>=1 or cor/=1;
            fails=fails+1;
         elseif ret==0 and cor==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

   if fails==0;
      success=0; bta=0;
      do while success==0 and fails<=30;

CMLset;
/* Alternative beta estimation */
__title=" ** A Beta Model, Numerical Solution (Alternative Estimation) **";

stvala=0.0|ones(k-1,1)|2.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvala=(bta')/2;
         print "retrying run";
      endif;

proc gradpa(theta,ds);
    local y,one,x,z,be,h,xb,zh,zh1,e,per1,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    zh1=zh+1;
    e=xb./(1+xb);
    per1=(-x.*e+x.*(e^2)).*zh1;

    gr1=psi(e.*zh1+(1-e).*zh1-1).*(x.*e.*zh1-(e^2).*zh1.*x+per1)-
        psi(e.*zh1-e).*(x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x)-
        psi((1-e).*zh1-1+e).*(per1+x.*e-(e^2).*x)+
        (x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x).*ln(y)+
        (per1+x.*e-(e^2).*x).*ln(1-y);
    gr2=psi(e.*zh1+(1-e).*zh1-1).*(e.*z.*zh+(1-e).*z.*zh)-
        (psi(e.*zh1-e).*xb.*z.*zh)./(1+xb)-
        psi((1-e).*zh1-1+e).*(1-e).*z.*zh+xb.*z.*zh.*ln(y)./(1+xb)+
        (1-e).*z.*zh.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglika(theta,ds);
    local y,one,x,z,be,h,xb,zh,e,v,p,q,llik,ds1;
		     
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    e=xb./(1+xb);
    v=e.*(1-e).*(1/(zh+1));
    p=((e^2).*(1-e))./v-e;
    q=(e.*(1-e)^2)./v-(1-e);
    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradpa;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400; 
    {ba,logl,g,vca,ret}=CMLprt(CML(data,0,&loglika,stvala));

/* check that parameters are reasonable */
   bta=(ba');

/* determine success or failure of the trial */
   cora=vca[1,1]/vca[1,1];
         if ret>=1 or cora/=1;
            fails=fails+1;
         elseif ret==0 and cora==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

/* linear heteroskedastic estimation */

   if fails==0;
      success=0; bth=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvalh=inv(x'*x)*x'*y|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=10;
         stvalh=(bth')/2;
         print "retrying run";
      endif;  

proc gradph(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikh(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ds[.,1];    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradph;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {blh,logl,g,vch,ret}=CMLprt(CML(data,0,&loglikh,stvalh));

/* check that Monte carlo parameters are reasonable */
bth=(blh');

/* determine success or failure of the trial */
   corh=vch[1,1]/vch[1,1];
       if ret>=1 or corh/=1;
          fails=fails+1;
          elseif ret==0 and corh==1;
          success=1; fails=0;
       endif;
    endo;
endif;

/* linear heteroskedastic estimation (with logit tranformation) */

   if fails==0;
      success=0; btl=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvall=inv(x'*x)*x'*(ln(y./(1-y)))|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvall=(btl')/2;
         print "retrying run";
      endif;  

proc gradpl(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ln(ds[.,1]./(1-ds[.,1]));
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikl(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ln(ds[.,1]./(1-ds[.,1]));    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradpl;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {bll,logl,g,vcl,ret}=CMLprt(CML(data,0,&loglikl,stvall));

/* check that Monte carlo parameters are reasonable */
btl=(bll');

/* determine success or failure of the trial */
   corl=vcl[1,1]/vcl[1,1];
      if ret>=1 or corl/=1;
         fails=fails+1;
      elseif ret==0 and corl==1;
         success=1; fails=0;
      endif;
   endo;
endif;

/* if trial has failed 30 times */

   if fails==30;
      print "run failed after 30 attempts";
      crash=crash+1;

      if ystatf==0;
         ystatf=(meanc(y)~stdc(y)~minc(y)~maxc(y));
      elseif ystatf/=0;
         ystatf=ystatf|(meanc(y)~stdc(y)~minc(y)~maxc(y));
      endif;

      if xstatf==0;
         xstatf=(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|maxc(x[.,2]));
      elseif xstatf/=0;
         xstatf=xstatf~(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|
                 maxc(x[.,2]));
      endif;
      fails=0;
      clear y;
   endif;

/* collect data for a successful trial */
   if success==1;

/* collect y stats for a successful trial */
      if trial==1;
         ystats=meanc(y)~stdc(y)~minc(y)~maxc(y);
         xstats=(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=b;
         sesum=sqrt(diag(vc));
         bsuma=ba;
         sesuma=sqrt(diag(vca));
         bsumh=blh;
         sesumh=sqrt(diag(vch));
         bsuml=bll;
         sesuml=sqrt(diag(vcl));
      endif;
      if trial>=2;
         ystats=ystats|(meanc(y)~stdc(y)~minc(y)~maxc(y));
         xstats=xstats~(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=bsum~b;
         sesum=sesum~(sqrt(diag(vc)));
         bsuma=bsuma~ba;
         sesuma=sesuma~sqrt(diag(vca));
         bsumh=bsumh~blh;
         sesumh=sesumh~sqrt(diag(vch));
         bsuml=bsuml~bll;
         sesuml=sesuml~sqrt(diag(vcl));
      endif;

/* vectors of first difference values for the j^th independent var */
         bns=b[1:cols(x),.];                @ nxk matrix @
         hns=b[cols(x)+1:rows(b),.];   @ nxk matrix @

/* calculate model mse */

         ahat=exp(x*bns);
	 bhat=exp(x*hns);
         yhatns=ahat./(ahat+bhat);
	 bnsmse=sumc((y-yhatns)^2)/n;
	 
/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alns=exp(xl'*bns);                          @ nxk-1 matrix @
         blns=exp(xl'*hns);
         ahns=exp(xh'*bns);
         bhns=exp(xh'*hns);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlns=alns./(alns+blns);                    @ nxk-1 matrix @
         meanhns=ahns./(ahns+bhns);
         varlns=(alns.*blns)./(((alns+blns)^2).*(alns+blns+1));
         varhns=(ahns.*bhns)./(((ahns+bhns)^2).*(ahns+bhns+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmns=meanhns-meanlns;
         fdvns=varhns-varlns;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmns=fdmns;  @ k-1X1 matrix for mean @
            totvns=fdvns;  @ k-1X1 matrix for variance @
	    totbmse=bnsmse;
         endif;

         if trial>=2 and ret==0; 
            totmns=totmns~fdmns;  @ k-1X(trial) matrix for mean @
            totvns=totvns~fdvns;  @ k-1X(trial) matrix for variance @
	    totbmse=totbmse|bnsmse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasim=(rndmn(b,vc,1000));  @ nxk matrix, where n=# of draws @
         be=betasim[.,1:cols(x)];                @ nxk matrix @
         h=betasim[.,cols(x)+1:cols(betasim)];   @ nxk matrix @

/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         al=exp(be*xl);                          @ nxk-1 matrix @
         bl=exp(h*xl);
         ah=exp(be*xh);
         bh=exp(h*xh);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlow=al./(al+bl);                    @ nxk-1 matrix @
         meanhigh=ah./(ah+bh);
         varlow=(al.*bl)./(((al+bl)^2).*(al+bl+1));
         varhigh=(ah.*bh)./(((ah+bh)^2).*(ah+bh+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdm=meanhigh-meanlow;
         fdv=varhigh-varlow;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totm=meanc(fdm);  @ k-1X1 matrix for mean @
            totv=meanc(fdv);  @ k-1X1 matrix for variance @
            totms=stdc(fdm);  @ k-1X1 matrix for mean @
            totvs=stdc(fdv);  @ k-1X1 matrix for variance @
         endif;

         if trial>=2 and ret==0; 
            totm=totm~meanc(fdm);  @ k-1X(trial) matrix for mean @
            totv=totv~meanc(fdv);  @ k-1X(trial) matrix for variance @
            totms=totms~stdc(fdm);  @ k-1X(trial) matrix for mean @
            totvs=totvs~stdc(fdv);  @ k-1X(trial) matrix for variance @
         endif;

/* get relevant stats for alternative method of beta estimation */
/* vectors of first difference values for the j^th independent var */
         bnsa=ba[1:cols(x),.];                @ nxk matrix @
         hnsa=ba[cols(x)+1:rows(ba),.];   @ nxk matrix @

/* calculate mse for alternative beta estimation */
  
         yhatnsa=exp(x*bnsa)./(1+exp(x*bnsa));
 	 bnsamse=sumc((y-yhatnsa)^2)/n;
	 
/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlnsa=(exp(xl'*bnsa)./(1+exp(xl'*bnsa)))';   @ nxk-1 matrix @
         meanhnsa=(exp(xh'*bnsa)./(1+exp(xh'*bnsa)))';
         dlsa=(1/(exp(xl'*hnsa)+1))';
         dhsa=(1/(exp(xh'*hnsa)+1))';
         varlnsa=meanlnsa.*(1-meanlnsa).*dlsa;
         varhnsa=meanhnsa.*(1-meanhnsa).*dhsa;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmnsa=meanhnsa-meanlnsa;
         fdvnsa=varhnsa-varlnsa;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmnsa=fdmnsa;  @ k-1X1 matrix for mean @
            totvnsa=fdvnsa;  @ k-1X1 matrix for variance @
	    totbamse=bnsamse;
         endif;

         if trial>=2 and ret==0; 
            totmnsa=totmnsa~fdmnsa;  @ k-1X(trial) matrix for mean @
            totvnsa=totvnsa~fdvnsa;  @ k-1X(trial) matrix for var @
	    totbamse=totbamse|bnsamse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasima=(rndmn(ba,vca,1000));  @ nxk matrix, where n=# of draws @
         bea=betasima[.,1:cols(x)];                @ nxk matrix @
         ha=betasima[.,cols(x)+1:cols(betasima)];   @ nxk matrix @

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanla=(exp(bea*xl)./(1+exp(bea*xl)));   @ nxk-1 matrix @
         meanha=(exp(bea*xh)./(1+exp(bea*xh)));
         dl=(1/(exp(ha*xl)+1));
         dh=(1/(exp(ha*xh)+1));
         varla=meanla.*(1-meanla).*dl;
         varha=meanha.*(1-meanha).*dh;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdma=meanha-meanla;
         fdva=varha-varla;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totma=meanc(fdma);  @ k-1X1 matrix for mean @
            totva=meanc(fdva);  @ k-1X1 matrix for variance @
            totmsa=stdc(fdma);
            totvsa=stdc(fdva);
         endif;

         if trial>=2 and ret==0; 
            totma=totma~meanc(fdma);  @ k-1X(trial) matrix for mean @
            totva=totva~meanc(fdva);  @ k-1X(trial) matrix for variance @
            totmsa=totmsa~stdc(fdma);
            totvsa=totvsa~stdc(fdva);
         endif;

/* set parameter vectors (linear heteroskeadastic) */
         bhh=blh[1:cols(x),1];
         hh=trimr(blh,rows(bhh),0);

/* calculate mse */

         yhatlh=x*bhh;
	 lhmse=sumc((y-yhatlh)^2)/n;
	 
/* get predicted means */
         elh=xl'*bhh;
         ehh=xh'*bhh;

/* get predicted variances */
         vlh=exp(xl'*hh);
         vhh=exp(xh'*hh);

/* produce first differences */
         if trial==1 and ret==0;
         totfdmh=ehh-elh;
         totfdvh=vhh-vlh;
	 totlhmse=lhmse;
         elseif trial>=2 and ret==0;
         totfdmh=totfdmh~(ehh-elh);
         totfdvh=totfdvh~(vhh-vlh);
	 totlhmse=totlhmse|lhmse;
         endif;

/* set parameter vectors (linear logit tranformation) */
         blo=bll[1:cols(x),1];
         hl=trimr(bll,rows(blo),0);
	 
/* calculate mse for linear het log trans */
  
         yhatlhlt=exp(x*blo)./(1+exp(x*blo));
	 lhltmse=sumc((y-yhatlhlt)^2)/n;

/* get predicted means */
         eol=exp(xl'*blo)./(1+exp(xl'*blo));
         eoh=exp(xh'*blo)./(1+exp(xh'*blo));

/* get predicted variances */
         vol=((exp(xl'*blo)./(1+exp(xl'*blo))^2)^2).*exp(xl'*hl);
         voh=((exp(xh'*blo)./(1+exp(xh'*blo))^2)^2).*exp(xh'*hl);

/* produce first differences */
         fdml=eoh-eol;
         fdvl=voh-vol;

/* produce first differences */
         if trial==1 and ret==0;
         totfdml=fdml;
         totfdvl=fdvl;
	 totthmse=lhltmse;
         elseif trial>=2 and ret==0;
         totfdml=totfdml~fdml;
         totfdvl=totfdvl~fdvl;
	 totthmse=totthmse|lhltmse;
         endif;

/* run ols */

{ vnam,m,bols,stb,vcols,stderr,sigma,cx,rsq,resid,dbw } = ols(0,y,x);

/* calculate mse */

      yhat=x*bols;
      olsmse=sumc((y-yhat)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsum=bols;
         selsum=stderr;
      endif;
      if trial>=2;
         blsum=blsum~bols;
         selsum=selsum~stderr;
      endif;

/* obtain the first differences */

olsl=xl'*bols;
olsh=xh'*bols;
fdols=olsh-olsl;

      if trial==1;
         fdolsum=fdols;
	 totomse=olsmse;
      endif;
      if trial>=2;
         fdolsum=fdolsum~fdols;
	 totomse=totomse|olsmse;
      endif;

/* run ols (with logit trans) */

yl=ln(y./(1-y));

{ vnam,m,bolsl,stb,vcols,stderrl,sigma,cx,rsq,resid,dbw } = ols(0,yl,x);

/* calculate mse */

      yhatlt=exp(x*bolsl)./(1+exp(x*bolsl));
      ltmse=sumc((y-yhatlt)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumo=bolsl;
         selsumo=stderrl;
      endif;
      if trial>=2;
         blsumo=blsumo~bolsl;
         selsumo=selsumo~stderrl;
      endif;

/* obtain the first differences */

olslo=exp(xl'*bolsl)./(1+exp(xl'*bolsl));
olsho=exp(xh'*bolsl)./(1+exp(xh'*bolsl));
fdolsl=olsho-olslo;

      if trial==1;
         fdolsuml=fdolsl;
	 tototmse=ltmse;
      endif;
      if trial>=2;
         fdolsuml=fdolsuml~fdolsl;
	 tototmse=tototmse|ltmse;
      endif;

/* run ols (with log trans) */

ly=ln(y+.01);

{ vnam,m,bolsll,stb,vcols,stderrll,sigma,cx,rsq,resid,dbw } = ols(0,ly,x);

/* calculate mse */

      yhatll=exp(x*bolsll)-0.01;
      ltmsel=sumc((y-yhatll)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumol=bolsll;
         selsumol=stderrll;
      endif;
      if trial>=2;
         blsumol=blsumol~bolsll;
         selsumol=selsumol~stderrll;
      endif;

/* obtain the first differences */

olslol=exp(xl'*bolsll)-0.01;
olshol=exp(xh'*bolsll)-0.01;
fdolsll=olshol-olslol;

      if trial==1;
         fdsumll=fdolsll;
	 totmsell=ltmsel;
      endif;
      if trial>=2;
         fdsumll=fdsumll~fdolsll;
	 totmsell=totmsell|ltmsel;
      endif;

/* reset fails and add a trial for a successful trial */
      fails=0;
      trial=trial+1; 
      clear y;

   endif;

endo;

/* check that parameters are accurate */
print "descriptive stats for standard beta parameters";
meanc(bsum')~stdc(bsum'); 
print " ";
/* check the mse for the std beta */
  
print "standard beta mse=";;meanc(totbmse);
print " ";
print "alternative beta mse=";;meanc(totbamse);
print " ";
print "ols mse=";;meanc(totomse);
print " ";
print "logit trans ols mse=";;meanc(tototmse);
print " ";
print "het ols mse=";;meanc(totlhmse);
print " ";
print "logit trans het ols mse=";;meanc(totthmse);
print " ";
print "log trans het ols mse=";;meanc(totmsell);
print " ";

/* check mean and std. dev of first differences */

print"first differences for mean and variance for std beta (sim)";
meanc(totm')~stdc(totm');
meanc(totv')~stdc(totv');
print " ";

print "first differences for mean and variance for std beta (nosim)";
meanc(totmns')~stdc(totmns');
meanc(totvns')~stdc(totvns');
print " ";

print "first differences for mean and variance for alt beta (sim)";
meanc(totma')~stdc(totma');
meanc(totva')~stdc(totva');
print " ";

print "first differences for mean and variance for alt beta (nosim)";
meanc(totmnsa')~stdc(totmnsa');
meanc(totvnsa')~stdc(totvnsa');
print " ";

/* check for mean and std. dev. of first difference for OLS */
print"first differences for mean ols";
meanc(fdolsum')~stdc(fdolsum');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (logit trans)";
meanc(fdolsuml')~stdc(fdolsuml');
print " ";

print "first differences for mean and variance for het normal";
meanc(totfdmh')~stdc(totfdmh');
meanc(totfdvh')~stdc(totfdvh');
print " ";

print "first differences for mean and variance for logit trans.";
meanc(totfdml')~stdc(totfdml');
meanc(totfdvl')~stdc(totfdvl');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (log trans)";
meanc(fdsumll')~stdc(fdsumll');
print " ";

/* efficiency of expected values */
/* compare efficiency of standard beta with alternative beta */

effnba=((totmnsa-fdmt)^2);
effdba=((totmns-fdmt)^2);
effalt=100*(sqrt(sumc(effnba')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. alternative beta=";;effalt;
print " ";

/* compare efficiency of standard beta with ols */

effno=((fdolsum-fdmt)^2);
effdb=((totmns-fdmt)^2);
effols=100*(sqrt(sumc(effno')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols=";;effols;
print " ";

/* compare efficiency of standard beta with ols */

effnl=((fdolsuml-fdmt)^2);
effdb=((totmns-fdmt)^2);
effolsl=100*(sqrt(sumc(effnl')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols (log trans)=";;effolsl;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbh=((totfdmh-fdmt)^2);
effdba=((totmns-fdmt)^2);
effhet=100*(sqrt(sumc(effnbh')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. het normal=";;effhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbll=((totfdml-fdmt)^2);
effdbal=((totmns-fdmt)^2);
efflog=100*(sqrt(sumc(effnbll')))./(sqrt(sumc(effdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;efflog;
print " ";

/* compare efficiency of standard beta with log trans normal */

effnbl=((fdsumll-fdmt)^2);
effdba=((totmns-fdmt)^2);
efflogl=100*(sqrt(sumc(effnbl')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. OLS (log trans)=";;efflogl;
print " ";

/* efficiency for variance term */
/* compare efficiency of standard beta with alternative beta */

veffnba=((totvnsa-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffalt=100*(sqrt(sumc(veffnba')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. alternative beta=";;veffalt;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbh=((totfdvh-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffhet=100*(sqrt(sumc(veffnbh')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. het normal=";;veffhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbll=((totfdvl-fdvt)^2);
veffdbal=((totvns-fdvt)^2);
vefflog=100*(sqrt(sumc(veffnbll')))./(sqrt(sumc(veffdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;vefflog;
print " ";

/* check standard beta parameters' overconfidence */

npoca=(bsum-meanc(bsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesum^2;
dpoc=sqrt(sumc(dpoca'));
ocp=100*(npoc./dpoc);
print "standard beta parms' overconfidence=";;ocp;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check alternative beta parameters' overconfidence */

npoca=(bsuma-meanc(bsuma'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuma^2;
dpoc=sqrt(sumc(dpoca'));
ocpa=100*(npoc./dpoc);
print "alternative beta parms' overconfidence=";;ocpa;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols parameters' overconfidence */

npoca=(blsum-meanc(blsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsum^2;
dpoc=sqrt(sumc(dpoca'));
oco=100*(npoc./dpoc);
print "ols parms' overconfidence=";;oco;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log. trans) parameters' overconfidence */

npoca=(blsumo-meanc(blsumo'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumo^2;
dpoc=sqrt(sumc(dpoca'));
ocol=100*(npoc./dpoc);
print "ols (log. trans) parms' overconfidence=";;ocol;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log trans) parameters' overconfidence */

npoca=(blsumol-meanc(blsumol'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumol^2;
dpoc=sqrt(sumc(dpoca'));
ocoll=100*(npoc./dpoc);
print "ols (log trans) parms' overconfidence=";;ocoll;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic normal parameters' overconfidence */

npoca=(bsumh-meanc(bsumh'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesumh^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic normal parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic logit transform parameters' overconfidence */

npoca=(bsuml-meanc(bsuml'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuml^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic log. trans. parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check overconfidence for the first differences (mean) */

print "standard beta's mean first difference overconfidence";
nmoca=(totm-meanc(totm'))^2;
nmoc=sqrt(sumc(nmoca'));
dmoca=totms^2;
dmoc=sqrt(sumc(dmoca'));
ocm=100*(nmoc./dmoc);
ocm;
print " ";

/* check overconfidence for the first differences (variance) */

print "standard beta's variance first difference overconfidence";
nvoca=(totv-meanc(totv'))^2;
nvoc=sqrt(sumc(nvoca'));
dvoca=totvs^2;
dvoc=sqrt(sumc(dvoca'));
ocv=100*(nvoc./dvoc);
ocv;
print " ";

/* check overconfidence for the first differences (mean) */

print "alternative beta's mean first difference overconfidence";
nmocaa=(totma-meanc(totma'))^2;
nmoca=sqrt(sumc(nmocaa'));
dmocaa=totmsa^2;
dmoca=sqrt(sumc(dmocaa'));
ocma=100*(nmoca./dmoca);
ocma;
print " ";

/* check overconfidence for the first differences (variance) */

print "alternative beta's variance first difference overconfidence";
nvocaa=(totva-meanc(totva'))^2;
nvoca=sqrt(sumc(nvocaa'));
dvocaa=totvsa^2;
dvoca=sqrt(sumc(dvocaa'));
ocva=100*(nvoca./dvoca);
ocva;
print " ";

/* check percentage of parameters within 95% conf int. */

print "check standard beta's percentage within 95% conf int";
cib=abs(bsum-meanc(bsum')).<(1.96*sesum);
sumc(cib')/(trial-1);
print " ";

print "check alternative beta's percentage within 95% conf int";
ciba=abs(bsuma-meanc(bsuma')).<(1.96*sesuma);
sumc(ciba')/(trial-1);
print " ";

print "check ols percentage within 95% conf int";
cio=abs(blsum-meanc(blsum')).<(1.96*selsum);
sumc(cio')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cit=abs(blsumo-meanc(blsumo')).<(1.96*selsumo);
sumc(cit')/(trial-1);
print " ";

print "check het ols percentage within 95% conf int";
cih=abs(bsumh-meanc(bsumh')).<(1.96*sesumh);
sumc(cih')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cil=abs(bsuml-meanc(bsuml')).<(1.96*sesuml);
sumc(cil')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cill=abs(blsumol-meanc(blsumol')).<(1.96*selsumol);
sumc(cill')/(trial-1);
print " ";

/* get information about percentage of significant parameters */
/* mean values */
print "check percentage of sig parms (standard beta)";
tstatbm=(abs(totm./totms)).>1.96;
tstatbv=(abs(totv./totvs)).>1.96;
sumc(tstatbm')/(trial-1);
sumc(tstatbv')/(trial-1);
print " ";

print "check percentage of sig parms (alternate beta)";
tstatba=(abs(bsuma./sesuma)).>1.96;
sumc(tstatba')/(trial-1);
print " ";

print "check percentage of sig parms (basic ols)";
tstatbo=(abs(blsum./selsum)).>1.96;
sumc(tstatbo')/(trial-1);
print " ";

print "check percentage of sig parms (logit trans ols)";
tstatbl=(abs(blsumo./selsumo)).>1.96;
sumc(tstatbl')/(trial-1);
print " ";

print "check percentage of sig parms (het ols)";
tstatboh=(abs(bsumh./sesumh)).>1.96;
sumc(tstatboh')/(trial-1);
print " ";

print "check percentage of sig parms (het ols - logit trans)";
tstatblh=(abs(bsuml./sesuml)).>1.96;
sumc(tstatblh')/(trial-1);
print " ";

print "check percentage of sig parms (ols -- log tranform)";
tstatbm=(abs(blsumol./selsumol)).>1.96;
sumc(tstatbm')/(trial-1);
print " ";


/* obtain information about y's */
print "descriptive stats for y's for successful runs";
meanc(ystats);
print "descriptive stats for y's for failed runs";
meanc(ystatf);

/* obtain information about x's */
print "descriptive stats for x's for successful runs";
meanc(xstats');
print "descriptive stats for x's for failed runs";
meanc(xstatf');

print "number of failed runs";
crash;
output off;